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Center for Turbulence Research-1987 activities (exlcuding the Summer Program) 


Background: 

The Center for Turbulence Research is a research consortium for fundamental study of 
turbulent flows. It is jointly operated by NASA- Ames Research Center and Stanford 
University. The Center became operational early in 1987 at a first year funding of 
$500,000. 

Administrative Matters: 

Currently, Parvdz Moin, John Kim, and William C. Reynolds are the executive officers of 
the Center. Moin is-the Director, Kim, the Ames coordinator, and Reynolds the program 
coordinator. The Center has a Steering Committee that meets regularly to act on the 
applications for Post-Doctoral Fellowships at the Center. The current members of the 
Committee are: 

D. Chapman (Stanford) 

S. Davis (Ames) 

J. Kim (Ames) 

P. Moin (Ames/Stanford) 

W. Reynolds (Ames/Stanford) 

M. Rubesin (Ames) 

The Center also has a high level Advisory Committee that meets annually and reviews 
the Center’s activities and accomplishments. The Committee members represent 
government, industry, and academia. The the names of the Advisory Committee 
members are attached. Their first meeting took place on March 20-21, 1988. The 
chairman of the Committee reports its findings to the Center for Turbulence Research 
officers, to the Ames Director, and to the Director of Aerophysics. 

An office assistant (Susan Hinton) was hired on February 1 to help with administrative 
matters. 

Technical Activities: 

The following provides a chronology of the events and individuals who participated in 
CTR. The individuals listed were appointed before December 1, 1987. 

Dr. Tsan Shih was appointed (on April 1) as a Post-doctoral Fellow. Dr. Shih’s 
speciality is turbulence modeling and specifically, one point closures. He has worked 
closely with Nagi Mansour of Ames utilizing the direct numerical simulation databases 
for testing turbulence models. He is currently working on the near wall behavior of 
turbulence models. Shih, on partial support from Air Force Office of Scientific Research, 
has been working with Parviz Moin on modeling of three-dimensional turbulent 
boundary layers. He has carried-out simulations and documented the statistical behavior 
of the channel flow subjected to imposition of pressure gradient or shear in the transverse 
direction. He is currently computing the Reynolds stress budget in this flow. 

Dr. Laurence Keefe (June 1987-Nov. 1987) has continued his work at Ames to measure 
the dimension of an attractor in low Reynolds number turbulent channel flow. Since 
December 1 he has been supported by AFOSR. The overall objective is to study the 
relevance and implications of the dynamical systems theory to open turbulent flows. 


Dr. Kevin Thompson (June 1987-present) (1/3 CTR support, 2/3 Ames’ Space Science 
Division) is studying turbulence in the early solar nebula. In this environment, turbulence 
is subjected to stratification, variable gravity, variable (Keplerian) rotation, and heat 
sources and heat loss by radiation. Thompson is currently writing and testing a code for 
simulation of compressible turbulence in three-dimensions. 

Dr. Julian Hunt (July 1987-August 1987)(Cambridge Univ.) spent an additional month 
beyond the summer program to continue work on space time correlations in 
homogeneous turbulence, on sel-similarity of two-point correlations in boundary layers, 
and on rapid distortion theory of near wall turbulence. 

Professor Paolo Orlandi (July 1987-Oct 1987) continued the work on development of 
an accurate finite-difference code for simulation of incompressible Navier-Stokes 
equations in generalized geometries during three of his four months’ stay. The code is 
now tested for several well known laminar cases, and should be completed during Prof. 
Orlandi’s visit in the summer of 88. 

Professor C. Benocci (Sept. 14, 1987-Oct. 10, 1987) is an assistant professor at the von 
Karman Institute in Brussels. His objective was to test phenomenological models for 
Lagrangian statistics using direct numerical simulation of forced isotropic turbulence. 
Although the generality of the results were hampered by the low Reynolds numbers in 
the simulations, some new results emerged. In particular, it was shown that contrary to 
the earlier assumptions, Lagrangian and Eulerian velocity auto-correlation are quite close 
to each other. Moreover, the autocorrelation curve does not exhibit the exponential shape 
in contrast to the solutions to the commonly used Lagrangian (Longevin) equation. It 
should be noted that Benocci used a 32x32x32 simulation at very low Reynolds numbers, 
and these results should be reexamined with higher resolution and Reynolds numbers. 

Dr. James Broadwel! (Oct. 87-Dec. 87) is a Sr. scientist at Cal-Tech. He is working 
closely with M. Rogers and R. Moser at Ames. The objective is to test Broadwell’s 
model for mixing. This model is based on organized structures in shear flows. This work 
has led M. Rogers t o dev elop a new code for time-developing mixing layer and jets and 
is a key element for CTR’s plans for strengthening the combustion and reacting flow 
research. 

Professor N. Etemadi (September 1987-August 1988) (1/3 CTR support, 2/3 Univ. of 
Illinois) is a mathematician from University of Illinois in Chicago. His field of speciality 
is Probability Theory. His appointment was based on the Center’s goal of generating new 
ideas in turbulence research. Prof. Etemadi has had no prior knowledge of fluid 
mechanics or turbulence. He has used the first half of his sabbatical in studying classical 
turbulence and turbulence terminology. It is expected that in the second half he will apply 
his expertise in probability theory to turbulence. 

Dr. George Karniadakis’s (Oct 1987-February 1988) speciality is in applica tion o f 
spectral element method to flows in complex geometries. During his tenure at CTR the 
objective was to apply the Spectral Element Method to direct simulation of turbulence on 
a wall with riblets. It has been determined experimentally that the flow over such a wall 
has reduced skin friction despite the increase wetted area. To this end, the laminar flow 
over a triangle was computed, for code validation and the results compared well with the 
experimental results. The 2-D spectral element code was modified to three dimensions 
and preliminary runs for the channel with riblets is initiated. It is expected that these 
calculations will be continued by Dr. Kamiadakis at MIT, and a graduate student will 
continue the work at CTR. 



I Drs. Thomas, R. Osborn and Dr. Hidekatsu Yamazaki (Oct. 16, 87-Nov. 16 87). The 

objective of these oceanographers from The Johns Hopkins University was to study 
encounter rates between planktonic particles in a turbulent fluid. The application is to the 
food web of oceanic plankton. Random walk was used to model planktonic motion 
relative to the turbulent flow environment which was simulated with a 64x64x64 
calculation of forced isotropic turbulence. As expected, for the cases with least energetic 
random walks, turbulence increased the contact rates between prey and predator, and 
with the most energetic random walks, effect of turbulence was negligible. The results to 
date are very preliminary. The work is being continued at Johns H opkin s, and it is hoped 
that its eventual dissemination will attract other oceanographers to CTR. 

Dr. Moon J. Lee (Nov. 1987-present) will continue his work on analysis of the effects of 
shear on turbulence. - A simulation of a shear-free turbulent boundary layer is planned to 
complement the homogeneous shear flow simulations. 

Dr. Jonathan Watmuff (Nov. 1987-present) will conduct an experimental investigation 
of turbulent boundary layers with adverse pressure gradient. The objective is to test the 
limits of Spalart’s assumptions in his direct numerical simulation code. Dr. Watmuff will 
bring a great deal of expertise in flow instrumentation to Ames’ Fluid Mechanics 
Laboratory. 

Graduate Student Research Assistants: 

P. Beaudan (Oct. 1987-) 25% RA. Spectral Element Method for complex geometries. 

* J. Neves (Oct. 1987-March 1988) 50% RA. Numerical simulation of an axial-flow over a 

cylinder. J. Neves will be supported by Office of Naval Research, starting April 1, 1988. 

S. Sorakayalpet (Oct. 1987-March 1988) 25% RA. Space-time characteristics of wall- 
pressure fluctuations. S. Sorakayalpet will continue on an ONR contract. 

K. Squires (Oct. 1987-Sept. 1988) 50% RA. Effect of particle loading on turbulence 
structure and modeling. 

M. Plesniak (Jan. 1988-June 1988) Experimental study of the effects of longitudinal 
curvature on mixing layers. 

M. Woronowicz (July 1987-Sept 1987) This study was sponsored by CTR during the 
summer of 1987 to determine whether particle simulation could be used to study low- 
speed turbulence. The method has been developed by Baganoff to simulate hypersonic 
flows. Several standard laminar flows were simulated with moderate success. This 
involved developing methods to measure such quantities as kinematic viscosity from the 
particle simulations. Due to some anomalies in la minar flow calculations, no attempt was 
made to simulate turbulent flows. Presently, CTR has no plan to continue supporting this 
effort. 
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Lagrangian Statistics in Homogeneous Isotropic 
Turbulence: Comparison Between Direct 
Numerical Simulation and Random Flight Models 

By C* Benocci 1 

1 Institut von Karman de Dynamique des Fluides 

Lagrangian techniques have found widespread application to the prediction and understanding 
of turbulent transport phenomena (Monin and Yaglom) and have yielded satisfactory results for 
different cases of shear flow problems (Durbin, 1983). However, it must be kept in mind that in 
most experiments what is really available are Eulerian statistics (Monin and Yaglom), and it is 
far from obvious how to extract from them the information relevant to the Lagrangian behavior of 
the flow; in consequence, Lagrangian models still include some hypothesis for which no adequate 
supporting evidence was up to now available. Direct numerical simulation of turbulence offers a new 
way to obtain Lagrangian statistics and so verify the validity of the current predictive models and 
the accuracy of their results. After the pioneering work by Riley (Riley and Patterson, 1974) in the 
70*8, some such results have just appeared in the literature (Lee et al, Yeung and Pope). The present 
contribution follows in part similar lines, but focuses upon two particle statistics and comparison 
with existing models (Durbin, 1980, Sawford and Hunt, 1985). 


l.The classical Lagrangian model 

In Lagrangian modeling, turbulent transport phenomena are simulated by the motion of “fluid 
particles.” The fundamental variables are the position of the particle X(i,t), where x is the initial 
point of the trajectory (release point) and its velocity V(x, t). Another quite useful quantity is the 
displacement vector Y(x, t) defined as: 


Y(i,0=X(®,t)-x 

Reynolds decomposition can obviously be applied to those quantities, giving, for example 


v(M) = v(M) + v'(M) 


from X, Y, and V the relevant statistical information could be extracted, the most important 
being the fluctuating displacement covariance tensor: 

the fluctuating velocity covariance tensor: 

B'y’ =V'(x,t)V'(x,t) 

and the Lagrangian integral time scale: 


T <L| = f R l^dt = 
Jo 


r VMM)V'(£,o) 

/ = 1 /2 5 1 /2 

Jo V' 2 (x,t) V\ 2 (x,o) 


where R) t L * is the velocity autocorrelation. 
Yaglom) show that: 


Dimensional and theoretical analysis (Monin and 


DjfW-t 2 


t « t (L > 
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dJ, L| ~ t t» T^' 

more refined formulations can be obtained, for simple flows only making assumptions about the 
shape of the autocorrelation function Rh* l K The latter is generally (Lee et al) assumed to have an 
exponential form: 

Using the form above, an analytical relationship can be found for the 2nd order moments on 
the displacement (the trace of the displacement covariance tensor) for the case of homogeneous 
stationary turbulence: 

■S'f 

where u'< is the fluctuating Eulerian velocity and must still be determined. To this end, it 
has generally been assumed that the Lagrangian autocorrelation be stronger than the Eulerian one, 
and a linear relationship of the type: 


has often been used (/? ~ .4 for the atmospheric surface layer). However, in the first numerical 
results (Riley and Patterson, 1974), the size and shape of the Lagrangian and Eulerian correlations 
were actually quite close. A similar trend is also shown by the most recent tests (Lee et al) for ho- 
mogeneous decaying turbulence. In view of this discrepancy between theory and available numerical 
results, investigation of the Lagrangian velocity autocorrelation and time scale is the most urgent 
task of direct numerical simulation. 

2. The random walk model of Lagrangian transport 

Simplified relationships of the type presented in the previous section are available only for some 
very simple cases of turbulent flow; for all the others, numerical solutions have to be sought. To 
this end, the turbulent Eulerian velocity field is taken as a random process and replaced by a known 
stochastic process of similar statistics. For a turbulent flow having a finite time scale and high 
turbulent Reynolds number, it is assumed that the turbulent acceleration is uncorrelated and can 
be simulated with a Gaussian random walk process. Under this hypothesis, the infinitesimal change 
of Lagrangian velocity dv over the time dt is given by the Lagrangian equation (Durbin, 1983): 



where <fW t is a step in a random walk process having a Gaussian pdf of mean 0 and variance dt. 
The Lagrangian equation automatically imposes the exponential form of the Lagrangian autocorre- 
lation. In practice a discrete process is used and the particle is advanced in time with an explicit 
technique: 

Xf = X*” 1 4- V” _1 At 
V? = V"- 1 + AV, 

where n is the current time step. The above technique has two main weak points: the first lies in 
the fact that the trajectory of each particle is generated by an independent random process and is, 
therefore, entirely uncorrelated with respect to other trajectories, while in reality particles moving in 
a correlated Eulerian velocity field have to be correlated with each other; the second (and related) is 
that the velocity field so generated does not respect continuity. The first difficulty is removed by two 
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particle models (Durbin, 1980, Sawford and Hunt, 1985) where the basic “event* is the correlated 
motion of one pair of particles. The two particles are moved with a law of the type: 

dXi = (aU' + 0U ")dt 
<tt 2 = (aU" + 0V')dt 

where a and are functions of the pair separation 5, defined as (Durbin, 1980): 

X L -Xg 
b “ 2 1 / 2 


and the Lagrangian length scale: __ 

L = (u' 2 ) 1 / 2 T w 

while the U' and U" velocities are respectively: 

TT' - V i +V 2 tT// _ Vx-V2 

U “ 2 1 / 2 ” 2 1 / 2 

Vi and V 2 being generated through two independent random processes obeying the Lagrangian 
equation. The issue of continuity can be addressed only by treating reverse dispersion problems (Saw- 
ford and Hunt, 1985) (i.e. finding the source which corresponds to a given marker’s distribution). 
This makes the approach unpractical for real prediction work. Comparison with direct numerical 
simulation should give an evaluation of the sise of the errors due to the unphysical features of the 
velocity fields. 

3. Direct simulation of Lagrangian motion 

In this initial phase, the study of Lagrangian statistics was limited to the case of isotropic homo- 
geneous stationary turbulence with no mean velocity. The corresponding Eulerian velocity field was 
generated with the Rogallo code (Lee and Reynolds, 1985); stationarity in absence of a mean shear 
gradient was achieved by adding a forcing term on the lowest wavenumber. Most of the simulations 
made were for a 32 3 mesh. The relevant quantities of the Eulerian flowfield are: 


v = 0.1 

u' 2 = 2.25 (deviation 5%) 


e = 8.5 


Therefore, the Kolmogorov microscale is: 

tj = (y ) 1/4 = 0.104 


and the Taylor microscale: 

A = (15^7^ j) 1 / 2 = 0.945 

corresponding to a turbulence Reynolds number R*e = 21.26. The separation between the two 
scales is, therefore, too small for the existence of a developed inertial subrange. The theoretical 
relationship: 

L< e > = A AReA 

lo 

gives a value of 1.07 for the Eulerian macroscale L <E * compared to about 1.25 for the simulation. 
The turbulent Reynolds number based upon the length scale lies therefore in the 25-30 range. It has 
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to be remarked that the effect of the forcing is to introduce a large scale perturbation which cannot 
be removed from the results: the Eulerian correlation does not, therefore, go down to 0 (Figure 1), 
and the overall concept of macroscale remains vague. Before passing to the analysis of the results, 
it must be remarked again that the Lagrangian equation model requires the turbulent acceleration 
to be uncorrelated i.e. the acceleration time T a scale to be much smaller than the velocity one, as: 

(ReL) -1/3 T (X '> 

This given for the present simulation: 


T 0 0.2T ,i > 

We are, therefore, at the very limit of the region of applicability of the model. The results of the 
present comparison cannot, therefore, be taken as final until confirmed by tests at higher Re’s. 

4 . Influence of the number of markers 

The 32 3 simulations were run with 4,096 particles; the number was chosen to keep the CPU time 
requirement to a manageable amount 5 sec per time step). As stochastic methods are relatively 
slow to converge (error size decreases with N“*/ 2 where N is the number of events), the first task was 
to verify whether the sample was good enough to give reliable results. To this end the full results 
were compared with the ones obtained by sampling one “event” (particle or pair of particles) out 
of two and out of four. Some relevant comparisons are presented in Figure 2 in the form of relative 
difference between the two results, i.e.: 

Result 1 — Result 2 
output ^R e8U it i Result 2)/2 

Figure 2 compares the predicted 2nd order moment of the displacement of single particles for 4,096 
and 2,048 event: maximum relative error is less than 10“ 3 , showing that the present distribution 
is adequate for single particle statistics. Figure 3 show the same comparison for the second order 
moments of the separation of a pair of particles (2,048 versus 1,024) events) and the increase in error 
is quite visible for small separation times; where the effect of the initial position of the marker is felt; 
for separation times higher than 2(c^ 4T^), the error is still low. As can be expected, much higher 
errors will be encountered when considering the separation velocity w = ds/dt, which is computed 
numerically from the separation; the separation velocity rms difference is shown in Figure 4. It can 
be seen from the figure that differences are of the order of 10“ 1 Therefore, results for the separation 
velocity are to be considered as qualitative only. Reducing again by half the number of samples, the 
errors increase by nearly one order of magnitude. It is somewhat surprising to find good convergence 
with so limited a number of events (most two particle simulations use 10 4 events (Sawford and Hunt, 
1985)). This is due to the fact that the sample is still generated from a correct velocity field. 

5. Influence of the initial distribution 

As previously mentioned, two tests were performed: the first case used a separation between 
markers of order of 2A and the second case used a separation of the order of A. Consistent with the 
effect of the number of markers on the results, differences in single-particle statistics between cases 
using different initial distributions were negligible. Differences in two-particle statistics, however, 
are affected by the initial distribution, with the difference becoming negligible after the particles 
motion becomes uncorrelated (Figure 5). In conclusion, the two cases can be used indifferently to 
discuss all the statistics with the exception of the ones related to separation. 

6. Single particle Lagrangian statistics 

As it was shown in a previous paragraph, the present simulations 1 and 2 can be taken to represent 
fully converged and equivalent solutions of the Lagrangian problem as far as one particle statistics 
are concerned. Therefore, all the pertinent data will be examined together in the present paragraph. 
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Based upon the average of six results (3 components * 2 cases), the Lagrangian time scale takes the 
value: 

T {L) =0.526 

with a 5% dispersion for individual data. The corresponding length scale L is given by: 

L = = 1.18 

quite close to the value of the Euler ian macroscale (see paragraph 4). The present results, there- 
fore, agree with the previous ones (Riley and Patterson, 1974, Lee et al, Yeung and Pope) in 
indicating that the Lagrangian and Eulerian correlations are quite close. Whether these results are 
of general value or are due to the low Reynolds number of the simulation is an open question which 
can only be answered by further tests. A remarkable result of the present tests is the discrepancy 
between the computed Lagrangian autocorrelation and the theoretical exponential shape (Figure 6). 
It can be shown that the numerical solution of the Lagrangian equation will yield the exponential 
form for the autocorrelation. Similar trends are made evident by a comparison of the displacement 
p.d.f. (Figure 19). As an obvious consequence, the results of the direct simulation do not fit well 
with the already discussed theoretical results: 


T ' 2 = 2u' 2 T ,i,2 ( — - 1 + 
v T< l > 

for the 2nd order moment of the displacement. As it can be observed in Figure 8, the numerical 
results diverge from the previous relationship as soon as the initial development Y' ~ t 2 is over, and 
become parallel to the theoretical curve only for t > 16T^ L) . By contrast, the numerical solution of 
the Lagrangian equation fits the theory perfectly (Figure 9). To conclude, it must be remarked again 
that the present calculations have been made at the very lowest edge of the Reynolds number region 
where the Lagrangian equation can be applied, and, therefore, the above discussed results are not 
enough to draw a conclusion upon the overall validity of the model. Surely it draws the attention 
over a quite peculiar behavior of the Lagrangian variables at low turbulent Reynolds numbers and 
puts in evidence the need for both tests at higher Re and further tests in the same range to develop 
a suitable understanding of such a regime. 

7. Two particle statistics 

As it was already remarked, the statistics of a pair of particles are potentially the ones most 
important for the understanding of the Lagrangian phenomenology. However, the results reported 
here can only be regarded as preliminary in nature. It has already been observed (paragraph 5) 
that the number of “events” (about 2,000) is still low and the convergence of the results not yet 
complete. Moreover, the prediction appears to be strongly affected by the nature of the forcing. This 
can be remarked by comparing the separation average velocity obtained from the direct prediction 
(Figures 10) with that corresponding to the two particle random flight model (Figures 11). It can 
clearly be seen that the patterns of the mean velocities are quite different; the ones predicted by 
the Lagrangian equation decrease in time (for large separation times) with a law which is smooth 
enough (considering the relatively limited number of markers ) and filling the expected law. 

On the contrary, the direct simulation yields a profile containing consecutive peaks, more or less in 
correspondence with the peaks of the Eulerian correlation (Figure 1). It appears that at large times 
when the two elements of the pair are far away and their relative motion is dominated by the energy 
carrying eddies, the influence of the forcing becomes dominant. Consequently, large differences 
exist between the separation predicted by the direct simulation and the Lagrangian equation. As 
an example, the 2nd order moments are shown for the direct simulation in Figure 12 and for the 
two particle model in Figure 13. The differences are evident, above all in the initial part where 
the behavior of the mean separation velocity is entirely different. Unfortunately the only available 
theoretical prediction of the separation moments is valid only for high turbulent Re and assumes, 
above all, the existence of an inertial subrange. Therefore, its predictions are, unsurprisingly, far 
from the numerical results (Figure 14). 


6 


C. Benocci 


8 . Conclusions 

The present report is only a first quite preliminary draft covering some of the results obtained 
during the stage of the author at the Center for Turbulence Research. Another part of the results has 
not yet been processed and will be discussed in the forthcoming second draft. Tentative conclusions 
are: 

• It is confirmed that Eulerian and Lagrangian macroscales are of the same sise, at least in the 
range of turbulent Reynolds number which can, at present, be tackled by direct simulation. 

• The predicted autocorrelation has no negligible differences with respect to the expected exponen- 
tial shape. 

• One particle statistics can be analysed in detail using the present results, but two particle statistics 
require more "events” (double at least) and a better understanding of the influence of forcing upon 
the large scales. 

• Turbulent Reynolds numbers at least 4 times higher than the present one should be reached to 
fulfill all the aims of the present investigation. 
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OVERVIEW 


A feasibility study was conducted during the Summer of 1987 to determine 
whether particle simulation could be used to study low-speed turbulence in a gas. 
Additional work was conducted during the Fall of the same year to attempt to resolve 
several questions that were raised. In trying to determine whether or not turbulence 
can be modelled in this manner, one must first establish that a suitably large Reynolds 
number can be reached in the simulation. This requires that we have the means to 
measure such quantities as kinematic viscosity and Reynolds number based strictly on 
the behavior of interactions of thousands of particles. This method, which was 
originally developed by Baganoff for studying hypersonic flows, was adapted to 
analyze a particular low-speed scenario using an IBM AT microcomputer. Thirty-two 
thousand particles were used, and computation was restricted to two dimensions. 

Stokes' First Problem was considered as a suitable scenario for simulation. In it, 
a laminar boundary layer develops in both space and time following a similarity scaling 
law, so the Reynolds number and the shear stress are time-dependent. A program 
was developed to obtain data on the velocity profile, temperature, and shear stress, 
and to find a means for computing a Reynolds number. A wide range of freestream 
velocities were tested (0.4 < M m < 2.0). The lower bound on the Mach number was 

dictated by achieving the highest velocity possible before compressible effects became 
significant. The importance of this will be explained later. The upper bound was 
chosen to see what would occur once compressibility became important. Despite the 
pronounced effects of thermal (random) motion of the particles at low fluid velocities, 
the simulated velocity profiles fit the theoretical curves nicely. Also, one would expect 
to find a constant coefficient of kinematic viscosity for incompressible flows, and in this 
analysis, every simulation gave the same numerical coefficient. 

Once the value of the kinematic viscosity was found, the Reynolds number (Re) 
was computed. Reynolds numbers, based on the boundary layer thickness, of 1 0 - 400 
were obtained in these simulations. (This quantity was computed at different time 
steps as well as for different freestream velocities.) For a flat plate, a generally 
accepted, similarly defined Re for the beginning of transition to turbulence is 1220. In a 



duct, transition begins at an Re of 2300, based on the diameter of the duct. Also, for 
Couette flow, Couette measured a critical Re of 1900, based on the distance between 
the plates (Hinze. Turbulence , second edition, p. 76). 

So, it appears that these simulations can achieve the Reynolds numbers 
associated with the onset of turbulence. Judging from the existing programs on the AT, 
one can estimate the time requirements needed when running such a program on the 
CRAY2 computer. While it takes roughly 30 seconds per time step using 32000 
particles on the AT, a code adapted to the CRAY2 could use one million particles 
undergoing time steps every 1.5 seconds. Such a program could easily simulate a 
three-dimensional flow. 


STOKES1 

It was decided to simulate a flow that could easily be verified by theory. If such a 
simulation was successful, one could generate more confidence in the method. An 
obvious, simple choice was to simulate the flow field of "Stokes' First Problem” (also 
known as the "Rayleigh Problem"). In the theoretical formulation (see Schlichting. 
Boundary Laver Theory , seventh edition, pp. 90-91), a semi-infinite fluid, bounded by a 
stationary wall, is initially at rest. At t = 0+, the wall is suddenly accelerated to a 
constant velocity of V w . In time, momentum is transferred to the fluid, and a 
time-dependent boundary layer develops. For this problem, the analytical solutions for 
the velocity profile and the shear stress are 

where q = -£= 

- v 4v * 

2 

t(x,t) « t w e where x w » 

One may define a Reynolds number with a length scale based on the distance 



2 



away from the wall where the velocity is 99% of the freestream velocity. The value of 
the similarity parameter ti for which erf ( tj ) = .99 is 1.8214. This gives us a 
time-dependent Reynolds number of 



In this program, it was easier to keep the wall stationary and have the fluid, 
initially having velocity everywhere, begin to slow down at the wall. During 

development of the program, = 0.6 was used to get the highest fluid velocity 
possible before compressible effects became significant. 


SCATTER 

One problem with using particle simulations for low speed flows is due to the fact 
that fluid motion usually consists of bulk motion of the flow superimposed on the 
random velocities of the individual particles. In a typical compressible flow, the bulk 
motion of the flow depends on the Mach number M, where 


Bulk motion a 

M« - 
a 

2 

where a * yRT 

Thermal motion 

a T a C 2 

2 

where C » 3RT 


So, a - C rms . Hence, as M increases, u / a - u / C rms increases, and thermal 
fluctuations become a smaller fraction of the total velocity. Also, since 


dispersion = y— 


where N = number of particles 
(normal distribution) 
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there is less statistical scatter with more particles in the region of interest and particle 
density increases. 


STOKES1 RESULTS 

In solving Stokes First Problem, one finds that the governing differential equation 
scales according to a similarity parameter x/Jav t . Because of this scaling, it was 
decided to save data of a particular simulation at time steps that were related to each 
other by a certain constant factor. Usually, data was saved for time steps from 32 to 
2048 that were related to each other by factors of two. 

The size of the computational box was 20 divisions high along the wall and 80 
divisions long normal to it. These divisions were used for calculating averaged values 
of physical quantities from the particle populations lying within each division. 
Thirty-two thousand particles were used in this two-dimensional system, because this 
number made good use of the allowable array lengths available with the large memory 
model on the AT. (Using the huge memory module would increase run times 
significantly.) As in the theoretical formulation, only the development of quantities 
normal to the wall were calculated. When data was to be obtained during a certain 
time step, the number of particles in each of the 80 columns was totaled, as well as the 
sums of each particles' three velocity components, their squares, and sums of 
x-velocity times y-velocity. Further, in computing velocity, density, temperature, 
pressure, and shear stress profiles normal to the wall, another averaging took place. 
In calculating a quantity, at each x position, corresponding to one column, the values of 
the two adjacent columns were added to it. The resulting value was then divided by 
the number density of those three columns. This data sampling occurred at a point in 
the program after the particles' positions had been updated but before the collision 
algorithm had been applied in order to get values before equilibrium occurred through 
collisions. 

Later, data from a run, which consisted of sets of data from various time steps, 


were called up, and velocity profiles were plotted. If these profiles were to follow the 
theoretical similarity scaling law, one could fix the time at an arbitrary value and vary 
the normal distance. Then, the profiles would collapse onto each other. It can be seen 
from the included figures that the data collapses on the similarity parameter very well. 
Figure 1 shows the similarity profile for 0.6 scaled to fit a time step t of 32. Figures 

2 and 3 show the same similarity profile scaled to t * 512 and 2048, respectively. 

At this point we may either pick a Reynolds number and then obtain a coefficient 
of kinematic viscosity, or we may find the kinematic viscosity and from that find the 
Reynolds number. The latter method was chosen, because one can easily fit the 
theoretical velocity profile to the scaled data and get a kinematic viscosity coefficient 
from the best-fitting similarity parameter. As expected, this coefficient remains constant 
regardless of the time step used to scale the data. Also, since the temperature profile 
is roughly constant (see figure 4), this coefficient remains constant over the entire 
range of Mach numbers tested. The best-fitting similarity parameter yielded a 
kinematic viscosity coefficient of .0562 squared divisions per time step. Because no 
mean free path was introduced in this model, this coefficient is not dependent on the 
density of particles in the system. After finding this coefficient, the Reynolds number Re 
and the wall frictional coefficient C f w were calculated. For the ranges of time steps and 

freestream velocities tested, Re ranged from 1 0.3 to 41 3 and Cf w ranged from 0.00985 
to 0.394. 

Once these quantities were found, the calculated dimensionless temperature and 
density profiles normal to the wall were plotted in order to determine whether or not 
incompressible assumptions were being violated. 

Then, the Mach number was increased to supersonic speeds. This was done in 
order to observe the effects of compressibility on the model and perhaps find a limit 
beyond which the velocity would no longer scale with the similarity parameter. Instead, 
it was discovered during this study that, although the temperature profile normal to the 
wall stays fairly constant at lower speeds, at higher speeds it builds up near the wall, 
and this build-up apparently scales with the same similarity parameter as the velocity! 
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Notice how the temperature data changes as the Mach number increases from 0.6 to 
2.0 at a fixed time scale of 2048 timesteps (figures 4 - 6). The temperature of the wall 
(T w ) was constrained to be the same as that of the initial, undisturbed fluid; but as the 

Mach number increases, there is increased collisional activity of the particles near the 
wall, so in that region, T increases. The result is the "bulge” in T which can clearly be 
seen in the supersonic data. Figures 6 through 8 show that for time step scales of 32 
through 2048, this bulge apparently scales with the similarity parameter x/ for 
fixed Mach number (M m = 2.0). Also, notice how the density profile in these figures 

changes with distance. As one might expect, the density for an ideal gas at constant 
pressure is lowest in the region where the temperature is highest. 

Finally, the theoretical dimensionless shear stress profile was compared to the 
shear stress data, as shown in figure 9. Although the function seems to fall within the 
scatter of the data points, there is too much scatter to tell how well it represents the 

data. We found this to be rather disappointing. At this time, it is not understood why 

/ 

there is so much variance in these data points. Attempts to resolve this question 
caused this report to be delayed. 


COMPUTATIONAL TIMES 

Computation on the AT took thirty seconds of user time per time step for 32000 
particles. Based on estimates done by Jeff McDonald when he adapted a similar code 
for the CRAY2 from the AT, running a modified version of this program on the CRAY2 
computer would achieve the same results for 10 6 particles in approximately 1.5 
seconds per time step. In order to do this, one would have to make use of vectorization 
on the CRAY2. This would allow about 2400 time steps per hour for a flow represented 
by over 31 times more particles than in this study. Moreover, such a flow could be 
three-dimensional. 

There is a tradeoff involved, however. If one wanted to re-create the model 
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studied this summer, one could take advantage of the fact that the dispersion, or 
statistical scatter, could be significantly decreased by using more particles in a 
two-dimensional box of identical dimensions. If 10 6 particles were used, the particle 
density would be over 30 times as great, and the dispersion would be over five times 
smaller. Alternately, one could keep the same particle density, roughly, and give the 
system a third dimension. For 10 6 particles, that third dimension would have 31 
divisions in length. So, a three-dimensional simulation has an accuracy disadvantage 
compared to its two-dimensional counterpart for a fixed number of particles. The 
resolution of this compromise would depend on the goal of the investigator. 


FREE SHEAR 

Also, a two-dimensional free shear flow was constructed, where two parallel 
particle streams of equal width, having the same velocity but travelling in opposite 
directions, were allowed to interact at some initial time. The simulation was displayed 
on a color monitor with one particle stream in red and the other in blue. As the system 
evolved in time, mixing of the streams was observed, with diffusional length growing 
with time. The overall velocity profiles, temperature T, and shear stress x , were 
computed from suitably averaged moments of velocity across the entire width of the 
system. As the speeds and particle densities were increased, there was less statistical 
scatter due to thermal motion and better visualization of the flow. Even though mixing 
was evident, the diffusional length grew with a square root dependence on time, just as 
it did in the Stokes First Problem model. This indicates that the mixing was laminar 
over the period observed (roughly 1000 time steps), which may not be surprising since 
this setup was not unlike that for Stokes First Problem. The difference with this setup 
was that the wall had been replaced by the boundary of a gas travelling in the opposite 
direction. 

At this time it is not clear whether or not such a flow would remain laminar or at 
some later time become turbulent. The length of time needed to observe significantly 
higher Reynolds numbers was prohibited by the time taken to run the program on the 
AT. 
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CONCLUSIONS 


So, it appears that particle simulation can model simple, low-speed, gaseous 
flows, and that useful physical quantities can be extracted from them. Unfortunately, at 
this time, the shear stress data has too much scatter in it to tell how well it is being 
predicted. In doing this exploratory study, the code was kept as simple as possible to 
get results quickly. Perhaps something important yet subtle has been omitted from the 
computation. Further work is being directed toward studying the shear stress. Also, it 
may be necessary to introduce a variable mean free path so that flow parameters such 
as the kinematic viscosity may be varied. An enhanced capacity to vary these 
parameters would increase the utility of this model to simulate different flows. 

At this time, it is not clear whether or not a particle simulation would be suitable 
for modelling turbulent flows. This would depend on the results obtained by making 
the changes discussed above. With a capability to increase Reynolds numbers 
efficiently, either through the use of a more powerful computer, or by an enhanced 
capability to vary the flow parameters, one could then simulate a flow which should 
become turbulent, resulting in a dramatic change in velocity and velocity moment 
profiles, and observe whether or not the particle model would exhibit such a transition. 
Such a test would better determine the suitability of this method for studying 
turbulence. 
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Select data_set no. 

/f 1 o *>py/U3o5 . 32 . s tksi 
/f I cppy/U3o5 . 64 . s tksi 
/f 1 opp«/U3©5 . 128 . s t ksl 
/f l oppy/U3o5 . 256 . s tksi 
/floppy/V3o5. 512, s tksi 
/f loppy /U3o5 . 1824 . s tksi 
f 1 o ppy/U3o5 . 2848 . s tksi 


Tin# scaled to t = 32 

Select 'sin. var. denon. 
4 * nu * t = 7.18 

nil = 8.8562 

M = 8.6888 

U® = 8.1888 
s 

Reynolds no. = 15.58 

C f = 8.262598 


Rgure 1 



S.l.ct data_s.t no. 

/r 1 oppy/U3o3 . 32 . s tksi 
/f loppy/U3o5 . 64. stksl 
/f 1 oppy/U3o5 . 128 . s tks 1 
/f loppy /V3o5. 236. stksl 
/f loppy /U3o3. 512 . s tksi 
/f loppy /W3o5. 1024. s tksi 
/f loppy /V3o5. 2940. s tksi 


Tim scaled to t = 512 

S.l.ct ' si n. var. d.noM. 
4 • nu * t = 115.00 
nu - 0.0562 
M_ = 0.6000 

U” = 0.1800 

33 

Reynolds no. = 61.99 

C r = 0.065650 


Rgure 2 
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Select d*ta_ set no. 

/£ I oppy/V3o5 . 32 . s tksl 
/f l oppy/U3o5 . 64. s tksl 
/floppy/U3o5. 128. stksl 
/f 1 oppy/U3o5 . 256 . s t ksl 
/f loppy/U3o5.512. stksl 
/floppy/U3o5. 1024. s tksl 
/f 1 oppy/U3o5 . 2048 .stksl 
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Select dat*_set no. 

/f ioppy/V3o5 . 32 . stksl 
/f l«ppy/U3o5 . 64. stksl 
/f 1 oppy/V3o5 . 128 . s tksl 
/l lopvy/Y3o5. 256. stksl 
/f 1 oppy/U3o5 . 512 . stksl 
/( \ oppy/U3o5 . 1024 . s tksl 
/f 1 oppy /U3o5 . 2048 . s tksl 
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Select data_set no. 

/f 1 op P y/U3o2 . 32 . stksl 
/f 1 oppy/V3o2 .64.s tksl 
! /f 1 oppy/V3o2 . 128 . s tksl 

/f 1 oppy/U3o2 . 256 . s tks 1 
/f loppy/V3o2 . 512. s tksl 
/floppy /U3o2. 1024. stksl 
/f 1 oppy/U3o2 . 2048 . s tksl 


x/jqp t < 4*nu*t ) 


Tine scaled to t = 2048 
= 1.5000 


Figure 5 



K/sqrtC4«nu*t> 


Select no. 

/floppy/Utwo32. stksl 
/f loppy /Ut*#o64. stksl 
/f loppy /U twol28 .stksl 
/f 1 oppy/V t«#o256 . s tksl 
/f I oppy/V t *©512 . s tksl 
/f loppy /Utwol024. stksl 
/f loppy /Utwo2048 .stksl 


Tine scaled to t = 2048 
M Q r 2.0000 


Figure 6 
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x/sqr t < 4«m**t ) 


Select no. Ti«e to t = 512 

/[ loppy/Utwo32. stksi M a = 2.0000 

/f 1 oppm/U two 64 . stksi 

/f loppy/V two!28 .stksi 

/t loppy /Utwo256. stksi 

/f loppy/V two512 . stksi 

/f Ioppy/U two 102 4 . s tksl 

/f 1 oppy/U t wo2048 .stksi 

Figure 7 



Select data_set no. 

/f 1 oppm/V two32 . stksi 
/f loppy /V two64. s tksl 
/f loppy/Utwol28 .stksi 
/f loppy /Vtwo256. stksi 
/f 1 oppy/U two512 .stksi 
/f l oppy /U two 1024 . s tksl 
/f 1 oppy/VJ t wo2048 . s tksl 


Time scaled to t = 32 

n m r 2.0000 
9 


Figure 8 
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x/sqpt(4«nu*t) 


Select data_se t no. 

/f I oppy/U fwo32 . stksl 
/f I oppy/U two£4. stksl 
/f I oppy/U two 128 .stksl 
/f 1 oppy/U t*o256 . stksl 
/f 1 oppy/U two5 12 .stksl 
/f 1 oppy/U t»ol924. stksl 
/f l oppy/U two2948 .stksl 


Tine scaled to t r 

= 2. MM 


Figure 9 
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Study of the Motion of Particles 
With a Random Walk in an Isotropic Turbulent Field 

Hidekatsu Yamazaki 
Thomas R. Osborn 
Chesapeake Bay Institute 
The Johns Hopkins University 


1. Introduction 

This report discusses some of the preliminary results obtained from our visit to Stanford 
Univ. during Oct.-Nov. 1987. The major objective of our research is to study encounter rates 
between planktonic particles in a turbulent fluid. The results of this study can be applied to 
the food web of oceanic plankton. The typical size and the motion of a single step of a 
zooplankter are comparable to the Kolomogorv scales for oceanic turbulence (where the rates of 
dissipation are 10* 1 ~ 10' 6 cm 2 /sec s ). Therefore we expect to observe some effects on 
planktonic motion induced by a turbulent field, and this has an influence on the contact rates 
between pray and predator. Although the intuitive argument is simple, the theoretical study is 
extremely complicated because of the spatial and temporal correlation of the turbulence field. 

Direct simulation of Navier-Stokes equation enables us to combine the effects of a 
turbulent velocity field with a random walk used to model planktonic motion relative to the 
water. We made use of pseudo-spectral code for turbulent simulations developed by NASA and 
the Stanford University group. An isotropic turbulent field was generated on 64 3 * * grid points 
using a forcing scheme. The forced simulation is desirable because we can study the variation 
in contact rate at a constant turbulent intensity. The velocity calculation started from a field 
that had reached a steady state. The kinematic viscosity was 0.1. The simulation showed that 
the Reynolds number was 20 and the rate of dissipation was 202. The time step interval was 
0.00046. We incorporated eight different sizes of random walks, in groups of 512 particles. 
The standard deviations of the random walks were the following, 

1) 0, namely no random walk, 

2 ) v, 

3 ) 

4) lOv 

5) 25 v 

6) 50v 

7) 100v 

8) 500v„ 

where v^ is Kolomogorv velocity scale (ev) 1 /* for a given rate of dissipation t and kinematic 
viscosity v. Before we pursued our major objective we examined the statistics of velocity 
vectors for each particle and the statistics of nearest neighbor distances. 

2. Particle contact rates 

We traced the eight groups of particles over 1200 steps (these calculation were a 

modification of K. Squires program, without whose help we could not have done this work in 

such a short time). The first 600 steps and last 600 steps were computed separately. We used 
the first 600 steps for the purpose of preliminary evaluation of contact rates between particles. 

Since we are interested in a difference in the contact rate with and without a turbulent field, 
we also simulated the motion of particles by a pure random walk with no turbulence. The 
particles were released from randomly distributed locations (by a uniform distribution) inside 
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the simulation domain (box). 


The calculation of contact was done in a second program after the program that moved 
the particles about. Particles that exited the box under consideration were moved back in 
through the opposite face so that the densities of predator or prey were not changed by 
particles escaping. When the distance between two particle becomes less than a threshold, the 
particles are considered to be in contact. We treated two cases for the prey particle after the 
contact. In the "without replacement" case, the particle (prey) was removed from the box. 
For the "with replacement" case, we would keep the particle "hidden" for a certain interval and 
then put it back into the system. The hiding of the particle was to avoid artificially large 
contact rates due to repeated contact between the same two particles. We present the result 
from the "without replacement" case in this report. It was found that the "with replacement" 
case had a source code problem and a corrected calculation is required. The with replacement 
calculation is the desirable result as it is the easiest to interpret. 


Since we have eight groups of particles many combinations of "prey-predator" conditions 
are possible. We calculated 28 cases of "prey-predator" contacts (Table 2.1). Predators with a 
large random walks contact all of prey before the 600th time step, so it is necessary to 
normalize the number of contacts. Let n be the number of prey contacted, n is less than 512 
for the small random walks, and N the total number of predators (512). The time step t, when 
the last prey contacted, can be used to average the time interval of searching for a prey. The 
normalized contact rates 7 is defined as follows, 


r\ 

7 = nT/(Nt) ; 


N 

rr r - 1 > ** »^*, r *■ *■ * ■ 


where T is the total time steps (600). The upper values of each cell is without turbulence and 
the lower value is with turbulence. The 6 cases for the smallest random walks reveal the 
result we expect, that turbulence should increase the contact rates between prey and predator. 
The rest of the values indicate that the turbulence field is an insignificant factor for the 
contact rates. While we are pleased with the preliminary result, it still needs to be put into a 
biological context (scaled with « and v) and then related to the analytical results in Rothschild 
and Osborn (attached and in press in the Journal of Plankton Research) 


Future work 

It is apparent that our results while encouraging are still primitive. There are both 
short term and long term aspects to the problem. In the short term we need to modify and 
correct our computer code that does the contact calculations. Some of this has been done 
already. We can then fully analyze the calculations from the 64 s model that were performed 
and saved this last fall. Those result can then be rationalized with the analytical calculations 
and form the basis for a first paper on the approach and the results. 

In the longer term, we need to perform some calculations at higher Reynolds number 
presumably on a 128 s grid. The planktonic motion needs to be expanded from a simple random 
walk to a series of flights that last over several time steps of the calculation. This change 
will increase the amount of calculation because another time scale is introduced (length of 
flight of the predator along the same path). As well, it will become important to incorporate 
behavioral responses to ingestion in the model. After eating, or when well fed, the predators 
do not hunt as hard; when starving, searching for food and the ability to escape predators is 
reduced. 

Essentially the short term projects are improvements of the technique and verification of 
the results. The long term projects are expanded development and utilization of the numerical 
calculation for studying the interaction of biological system with the turbulent fluid. 
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Table jil Contact rates between particles. The heading in each column and row refers to the 
random walk group. The upper number in each pair of data values is without random walk and 
the lower number is for random walk plus the turbulence. 


prey\predator 2 

3 

4 

5 

6 

7 

8 

1 

0.09 

0.17 

0.33 

0.77 

0.98 

1.39 

5.08 


0.25 

0.34 

0.46 

0.78 

0.99 

2.23 

6.67 

2 

X 

0.16 

0.32 

0.73 

0.99 

2.21 

6.12 


X 

0.31 

0.47 

0.79 

0.99 

1.57 

5.22 

3 

x - 

X 

0.32 

0.76 

0.98 

1.44 

5.88 


X 

X 

0.46 

0.78 

0.98 

2.46 

6.74 

4 

X 

X 

X 

0.76 

0.98 

1.73 

6.32 


X 

X 

X 

0.79 

0.98 

1.83 

5.45 

5 

X 

X 

X 

X 

0.99 

2.13 

4.48 


X 

X 

X 

X 

0.99 

1.97 

8.11 

6 

X 

X 

X 

X 

X 

2.33 

4.29 


X 

X 

X 

X 

X 

2.05 

5.41 

7 

X 

X 

X 

X 

X 

X 

5.17 


X 

X 

X 

X 

X 

X 

5.61 


3. Statistics of particles 


3.1 Velocity vectors of particles 

Velocity components generated from the turbulent simulation for each particle compose a 
set of time series. A velocity component at a particular time step correlates with the previous 
time step. Following the conventional notion of Lagrangian auto-correlation, the velocity 
correlation decreases exponentially as the time lag increases. Since we add a random walk to 
the particle motion at each time step, the conventional arguments are no longer valid. We 
examined the averages of the velocity time series of each particle. As the size of the random 
walk becomes larger, at each time step the particle is forced away from a position where the 
particle would be located by the turbulent field and the time series of velocity becomes 
uncorrelated with its own past. Therefore we expect the average values of the velocity 
components for large random walk should follow the distribution of the mean from white noise. 
On the other hand the average velocity components with no random walk must follow the 
distribution law of a correlated random variable. Let Uy, Vy and Wy be velocity components 
of i-th particle at j-th time step. These values comprise the velocity vector a and the 
magnitude of a is 


We examined the averages of velocity components u,v,w and velocity vector a over 200 time 
steps, namely the following values were calculated 
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and 


u r E Uy / N 
Vi* E Vy / N 
w,- E Wy / N 
A,« E ay / N. 

Table 3.1 -3.4 show the statistics of each average value. Since the random walk forces a 
particle to move in an uncorrelated fashion, the particle with the large steps tends to lose 
correlation from one time step to the next. The mean of U, V and W approach zero for larger 
random walks and the range of the statistics becomes smaller. A random sampling from the 
turbulence field, so that each time series of u, v and w were uncorrelated, showed a good 
agreement with case 8). Some of particles without random walk are trapped in eddies and some 
are scattered around. Therefore the averages of U, V and W showed the large standard 
deviation for a short realization period in the present case (200 steps). A careful inspection of 
the standard deviation indicates some peculiar interaction of turbulence field with random walk 
for cases 2) and 3). All of the statistics are monotonici with increasing random walk step size 
except these two cases. 

The random walk was generated from a normal distribution for all of our simulations. To 
ensure that no programing errors existed in the codes, a random walk based on the uniform 
distribution ( from the Cray intrinsic function) was also applied to each step and the result 
showed the same tendency as the normal distribution. 

Table 3.1 Statistics for Uj 


mean 

s.d 

min 

max 

1) -0.2405701 

4.4174653 

-11.1525369 

11.2729212 

2) 0.1061830 

4.6179897 

-12.4469889 

10.8448361 

3) -0.1063024 

4.4741023 

- 12.7496693 

10.6039775 

4) 0.1637198 

4.3686518 

-10.7305665 

10.0886646 

5) 0.0739651 

3.9530263 

-10.3270934 

11.1192905 

6) -0.0484558 

3.0176252 

-9.6488662 

8.6306654 

7) 0.0049584 

1.9290495 

-6.3919247 

5.7244361 

8) -0.0057150 

0.4719053 

-1.4193553 

1.2998666 


Table 3.2 Statistics for V t 


mean 

s.d 

min 

max 

1) -0.5182599 

4.1638463 

-12.6059994 

11.6511822 

2) -0.2592903 

5.0128338 

-14.2338938 

11.9635239 

3) 0.5663336 

4.6510084 

-10.7875769 

12.6121930 

4) 0.7882887 

4.5643164 

-9.9557783 

14.6300351 

5) -0.0366833 

4.1046428 

-9.7325182 

10.4931608 

6) 0.1122861 

2.9682378 

-9.5482282 

8.0167097 

7) 0.0345716 

1.9081249 

-6.4678161 

7.8802499 

8) -0.0079068 

0.4806509 

-1.8767631 

1.5038711 
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Table 3.3 Statistics for W s 


mean 

s.d 



1) -0.8477828 

4.7996655 

-11.5839232 

12.2917635 

2) -0.1974391 

4.6724461 

-14.8436250 

10.1079898 

3) -0.4030535 

4.9459295 

-12.9826462 

12.6153774 

4) -0.1622836 

4.7454426 

-14.2204286 

12.1482159 

5) 0.7569673 

4.2470180 

-13.7572350 

9.3310843 

6) 0.7373272 

3.0956171 

-8.6091972 

9.6507205 

7) 0.0727480 

1.8895971 

-6.1635830 

6.2546522 

8) -0.0084084 

0.5023862 

-1.5490949 

1.4636020 


Table 3.4 Statistics for A; 


mean 

$,d 

min 

max 

1) 7.8896638 

2.8654964 

0.9208719 

14.9406667 

2) 8.5395250 

2.9846524 

1.4578940 

16.9043228 

3) 8.5713267 

2.7004234 

2.3444629 

15.5812644 

4) 8.3468284 

2.8001345 

1.2910205 

16.5322217 

5) 8.3716425 

2.3055518 

2.8363484 

16.0422414 

6) 8.0586156 

1.6751876 

3.5614937 

12.8514342 

7) 7.9140421 

1.0677454 

4.8263259 

12.2567512 

8) 8.1199618 

0.3039232 

7.2648240 

9.1066186 


3.2 Statistics of nearest neighbor distance 

The distance between particles is important in our study. We examined the nearest 
neighbor distance among groups ( each has 512 points). A randomly distributed finite number 
of points in a finite range of space can be described by the theory of "point process". In 
general, the locations of particles in a field follows a Poisson process. The main objective of 
this section is to verify that the general theorem of nearest neighbor distance holds for 
particles in a turbulence field with/without random walks. The probability density function of 
the nearest neighbor distance r for the three dimensional case has the following form (Stoyan, 
Kendall and Mecke, 1987, p49). 

f(r) «= 3Akr 2 exp[-Akr s ] 

where k - 4 jt/ 3 and A is the intensity of poisson process. For a known number of points n and 
finite volume a, the value can be expressed 

A = n/a 

We used the box size of (2x) 3 for the simulation and each group contains 512 points. Thus 
A = 51 2/(2 jt) 3 

The expected value of r can be expressed 
E[r] - Jrf(r)dr 
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= (Ak)- 1 /*r(4/3) 

- 0.4354 

If the particles follow the general Poisson process, the average value of nearest neighbor 
distance should agree with the above value and the average value should not depend on 
whether it is with a random walk or a pure turbulent motion. We started with the particles 
randomly distributed (uniform distribution) in space for each case. The statistics of nearest 
neighbor were calculated at every 10 th step. Not only the average values (solid line) but the 
standard deviation (broken line), the minimum (triangle) and the maximum (circle) did not 
change from the beginning to the end for all cases (Figure 1). This confirms that particles as 
a point process among the same group follow the stationary Poisson process. 

4. Auto-Regressive models 

It is desirable to find a usable stochastic model for the turbulence velocity field to 
implement various types of motion of an organism in difference environments. The successful 
model would enable us to reduce the amount of computation and to consider larger ranges of 
spatial scales. We are not looking for a "perfect" turbulence field but rather one that mimics 
the statistics as far as particle encounters are concerned. This result would make simple 
modelling schemes available to a large number of biological modelers. 

A conventional approach to model the turbulence velocity field is Langevin’s equation 
which has a form of first order stochastic differential equation. A discrete form of the 
equation can be expressed as a first order discrete AR process. Many attempts have been made 
to simulate turbulence, particularly for turbulent diffusion. Although many authors claim the 
simulations give reasonably good agreement with the observations, they have not tested the 
model against the real turbulent velocity data. We attempted to fit many different orders of 
AR processes to the velocity data. Two criteria were used to identify the best fit to a given 
data. AIC is the most common order criterion (Akaike, 1971) but seems to give higher order 
than the true process (Priestley, 1981). Another criterion used is LIL (Hannan and Quinn, 
1979) and seems to give better results. 

Table 4.1 (AIC) and 4.2 (LIL) show the number of best order fitted to velocity component 
u for the eight different cases. Case 1 (no random walk) appeared higher AR models to fit the 
data than what we initially expected. However, this can be explained from the following 
reason. Suppose the turbulence velocity field is in a class of second order stochastic equation, 
the discrete data from such kind of continuous process tend to appear as Autoregressive/moving 
average process of order (2,1) ( ARMA(2,1) ). Since we are trying to fit AR model to ARM A 
process we should expect a high order of AR model for the estimate. The result suggests that 
Langevin’s equation (discrete form, so AR(1)) does not describe the velocity field well. The 
model may fit reasonably to the observed diffusion data but this may be an artifact of a 
statistical nature. 

As far as the finding the best model for the turbulence velocity field concerns that we 
should consider fitting ARMA models but we have not a chance to perform a test. 
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Table 4.1 The number of best AR order by AIC. The data is velocity component u. 
case 1 2 3 4 5 6 7 8 


1th order # of aic-= 

0 

0 

104 

216 

241 

231 

262 

325 

2th order # of aic= 

0 

4 

37 

57 

59 

63 

80 

85 

3th order # of aic- 

0 

8 

20 

32 

38 

44 

47 

37 

4th order # of aic= 

3 

4 

32 

40 

35 

34 

27 

10 

5th order # of aic* 

1 

6 

29 

31 

24 

38 

21 

11 

6th order # of aic= 

3 

4 

34 

19 

19 

20 

22 

12 

7th order # of aic- 

3 

8 

29 

20 

16 

16 

8 

7 

8th order # of aic= 

16 

43 

39 

19 

11 

17 

13 

7 

9th order # of aic- 

20 

37 

22 

10 

11 

13 

9 

4 

10th order # of aic- 

30 

68 

29 

16 

10 

12 

6 

3 

11th order # of aic= 

49 

71 

26 

11 

7 

8 

5 

2 

12th order # of aic= 

66 

78 

22 

11 

8 

5 

3 

2 

13th order # of aic* 

81 

60 

27 

8 

14 

5 

2 

3 

14th order # of aic= 

95 

50 

32 

9 

10 

2 

1 

0 

1 5th order # of aic= 

29 

24 

9 

7 

9 

2 

5 

3 

16th order # of aic= 

15 

10 

4 

2 

0 

1 

0 

0 

17th order # of aic= 

10 

2 

4 

1 

0 

0 

0 

0 

18th order # of aic= 

15 

12 

2 

0 

0 

1 

0 

1 

19th order # of aic= 

9 

4 

1 

0 

0 

0 

0 

0 

20th order # of aic= 

11 

4 

0 

1 

0 

0 

0 

0 

21th order # of aic= 

12 

6 

3 

0 

0 

0 

1 

0 

22th order # of aic= 

12 

1 

2 

0 

0 

0 

0 

0 

23th order # of aic= 

9 

5 

2 

1 

0 

0 

0 

0 

24th order # of aic= 

11 

0 

1 

0 

0 

0 

0 

0 

25th order # of aic= 

6 

2 

1 

0 

0 

0 

0 

0 


Table 4.2 The number of best AR order by LIL. The data is velocity component u. 
case 1 2 3 _4 5 6 7 g 


1th order # of lil= 

0 

1 

229 

366 

390 

384 

395 

434 

2th order # of lil= 

0 

11 

58 

51 

58 

64 

72 

56 

3th order # of lil= 

2 

14 

22 

27 

21 

24 

25 

13 

4th order # of lil= 

8 

5 

36 

21 

18 

10 

8 

3 

5th order # of lil= 

1 

9 

38 

22 

10 

11 

4 

3 

6th order # of lil= 

4 

19 

25 

6 

5 

4 

4 

1 

7th order # of lil= 

15 

37 

25 

5 

5 

3 

0 

1 

8th order # of lil= 

38 

94 

31 

7 

1 

7 

1 

1 

9th order # of lil- 

41 

49 

9 

1 

2 

3 

2 

0 

10th order # of lil- 

47 

77 

13 

1 

1 

0 

0 

0 

11th order # of lil- 

68 

59 

9 

2 

0 

2 

0 

0 

12th order # of lifr 

68 

61 

7 

0 

0 

0 

0 

0 

13th order # of lil- 

64 

28 

4 

0 

0 

0 

0 

0 

14th order # of lil- 

56 

22 

3 

2 

0 

0 

0 

0 

15th order # of lil= 

33 

12 

3 

1 

1 

0 

0 

0 

16th order # of lil- 

20 

6 

0 

0 

0 

0 

0 

0 

17th order # of lil= 

8 

1 

0 

0 

0 

0 

0 

0 

18th order # of lil= 

10 

3 

0 

0 

0 

0 

1 

0 

19th order # of lil= 

7 

2 

0 

0 

0 

0 

0 

0 

20th order # of lil- 

6 

0 

0 

0 

0 

0 

0 

0 

21th order # of lil= 

4 

1 

0 

0 

0 

0 

0 

0 

22th order # of lil= 

5 

0 

0 

0 

0 

0 

0 

0 

23th order # of lil- 

4 

1 

0 

0 

0 

0 

0 

0 

24th order # of lil= 

0 

0 

0 

0 

0 

0 

0 

0 

25th order # of lil- 

2 

0 

0 

0 

0 

0 

0 

0 
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Review of the Literature on Lagrangian modelling 

There exists an abundant volume of literature about the Lagrangian method in the field of 
fluid dynamics. Much is based on theoretical works. A limited number of laboratory 
experiments have been conducted by actually tracing a particle in space. Jones et al. (1967) 
made an early experiment to study the three-dimensional motion of particles. Sullivan (1971) 
and Snyder and Lumley (1971) made similar experiments, while more recently Well and Stock 
(1983) and Sato and Yamamoto (1987) used an optical device to trace a particle in three 
dimension. 

Numerical simulations have been applied in turbulence studies. The Lagrangian statistics 
are usually studied with "random walk" or Markovian processes. In his early paper on 
Lagrangian variables, Obukhov (1959) suggested that the evolution of a diffusion process can be 
described by a Markov process and obeyed the Fokker-Planck equation. The Langevin equation 
(Krasnoff and Peskin, 1971; Smith, 1968; Gifford, 1982) has been adapted to turbulent diffusion 
problems. Numerical simulations based on the Langevin equation were conducted by Reid (1978) 
and Hanna (1979). Ratm and Svensson (1986) applied the Langevin equation to oceanic diffusion 
within the Ekman layer, and obtained agreements with field observations of dispersion. 
Although simulation and observations show reasonable agreement, the Lagrangian 
auto-correlation function of velocity was in a form of an exponential function, and this 
constrains certain classes of turbulence from following the Langevin equation, particularly at 
low Reynolds number. 

Random-walk theory has been used in this study of animal motions, as early as Patlak 
(1953a,b). Fraenekl and Gunn (1961) classified undirected response (kineses) into orthokinesis 
(modulation of velocity) and klinokinesis (modulation of turning rates). Rohlf and Davenport 
(1969) discussed the effect of a stimulus field on animal movement pattern. Doucet and Drost 
(1985) categorized the movement in depth. Terrestrial animal (insect) movements have been 
conjugated with the population dynamics (Jones, 1977). Kareiva and Shigesada (1983) extended 
the model to the first order Markov process. Siniff and Jessen (1969) examined the movement 
of red foxes and showed that the mean speeds associated with each path follow the gamma 
distribution. The distribution of turning angles appeared to be normally distributed. Equivalent 
aquatic observation have not been reported in the literature. Sakai (1973) is the only study, at 
least known to the authors, on a simulation of fish grouping. The model is developed from the 
dynamic equation of motion with "forces" including forward thrust, mutual interaction and 
arrayal forces, but the random forcing term is rather deterministic. The simulated grouping 
pattern shows amoebic, toroid and rectilinear movements. 

Purcell (1977) and Zaret (1980) discuss the physical environment of planktonic organisms. 
In general the organisms live in highly viscous conditions at low Reynolds number. The 
population, however, is subject to turbulent diffusion and stirring. Accumulating evidence 
shows that planktonic food webs are strongly related to turbulent mixing (Sonntag and Parsons, 
1979; Gallegos and Platt, 1982; McGowan and Hayward, 1978). Sverdrup (1953) defined the 
concept of critical depth for the relationship between mixing layers and phytoplankton growth. 
Numerical model studies were conducted by applying Eulerian type of equations (e.g. Tett, 
1981). Recent review papers by Denman and Powell (1984) and Tell and Edwards (1984) discuss 
the effects of physical processes on planktonic organisms. Denman and Powell (1984) note that 
the relation of physical process and planktonic ecosystem should be studied by matching both 
biological and physical scales. By this account not only measurements of a planktonic 
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ecosystem but the modeling should be performed at the same order of scale. Eulerian types of 
models are based on averaged values over a certain physical scale, which is usually much larger 
than the size of the plankton. Tett and Edwards (19 ) mention that "Conversely, selection 
may operate to ensure that ’the plankton’ tends always to exploit the available resources, and 
in this case continuum model (Eulerian model) might be adequate for predicting total biomass 
and production..”. But they point out the Eulerian models are not suitable for studying the 
transient state and individual species of plankton. Then a Lagrangian model should be adopted. 
Lewis et al. (1984) apply an Eulerian model to study the mixing and photadaptation of 
phytoplankton, in which they claim that the Eulerian equations can be transformed into a 
Lagrangian coordinate system (they refer to Hill, 1976). This is "mathematically" true, but it is 
not tractable to obtain the Lagrangian statistics. 

Since the Lagrangian models track each individual particle (or plankter), extensive 
computational effort is required. But recent rapid advancements in computer technology 
decreased the cost of computer operation drastically. A numerical simulation model of 
Ledbetter (1979) is the earlier attempt of the Lagrangian model to study plankton patchiness in 
the Langmuir circulations. Similar numerical experiment (Watanabe and Harashima, 1986) shows 
good agreement of the Lagrangian model with the aggregating pattern with the Langmuir, 
circulation observed in a laboratory tank. Falkowski and Wirick (1981) use simple "random 
walk" to trace phytoplankton motion in the mixing layer. Woods and Onken (1982) use "random 
walk" plus sinusoidally modulated motion to study the effect of diurnal variation of mixing 
layers on the phytoplankton production. They follow a single plankter and run the same 
simulation 100 times to get an ensemble average of the plankton motion. Therefore the 
plankters are statistically and physically independent. Wolf and Woods (1987) extend Woods and 
Onken’s model (1982) to investigate further the effects of diurnal and seasonal variation of 
mixed layer depth on phytoplankton production. They suggest that at least 20 particles in each 
one-meter depth interval over the entire mixed layer, 5000 particles for their case, are 
necessary to ensure a statistically significant result. They show that physiological adaptation 
of plankton depends on the history of their individual trajectory. A hybrid Lagrangian model 
with a two-layer mixed layer model (Eulerian) examines the sinking characteristics of 
planktonic organism (Lande and Wood, 1987). 
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